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Abstract 



A numerically effective procedure for determining weakly reversible 
chemical reaction networks that are linearly conjugate to a known re- 
action network is proposed in this paper. The method is based on 
translating the structural and algebraic characteristics of weak re- 
versibility to logical statements and solving the obtained set of linear 
(in)equalities in the framework of mixed integer linear programming. 
The unknowns in the problem are the reaction rate coefficients and the 
parameters of the linear conjugacy transformation. The efficacy of the 
approach is shown through numerical examples. 

Keywords: chemical kinetics; stability theory; weak reversibility; linear 

programming; dynamical equivalence 
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1 Introduction 

A chemical reaction network is given by sets of reactants reacting at pre- 
scribed rates to form sets of products. The mathematical study of such 
networks has been applied recently to such fields as industrial chemistry, 



systems biology, gene regulation, among others [10,19,28 . There has also 
been significant theoretical work in the literature on such questions as per- 
sistence [lj|3] , multistability [8pp6] , monotonicity 4J[5] , the global attractor 



conjecture for complex balanced systems [T|[2|[7|[T3 , and conjugacy of reac 



tion networks 11,21 



One problem which has attracted recent attention has been that of de- 
termining when two reaction networks can exhibit the same qualitative dy- 
namics despite disparate network structure. It has been long known that 
two networks can given rise to the same governing set of differential equa- 



tions under the assumption of mass-action kinetics [20,22 . In [111 and [29], 
the authors complete the question of what network structures can given rise 
to a set of governing differential equations. This work is extended in f21j to 
networks which do not necessarily have the same set of differential equations 
but rather have trajectories related by a non-trivial linear transformation. 

The problem of algorithmically determining when a network has the 
same governing dynamics as another network satisfying specified proper- 
ties was first addressed in [30 j where the author presented a mixed- integer 
linear programming (MILP) algorithm capable of determining sparse and 
dense realizations, i.e. networks with the fewest or greatest number of re- 
actions capable of generating the given dynamics. This line of research was 
continued in |31| in which the authors extended the algorithm to complex 



and detailed balanced networks, and in [33| in which the authors addressed 
the problem of finding dense weakly reversible realizations. 

In this paper we show how the problem of determining weakly reversible 
realizations presented in [33] can be reformulated as a linear constraint 
within the established MILP framework. This reformulation significantly 
eases the computational cost associated with the problem which had previ- 
ously been solved through successive applications of an optimization algo- 
rithm and checking for weak reversibility with one of Kosaraju's, Tarjan's or 
Gabow's algorithm [6 23 . We also extend the algorithm to encompass the 



linearly conjugate networks introduced in 21 . We show how the algorithm 
can be used to effortlessly reproduce results from the literature and easily 
handle large-scale networks not yet considered. 

2 Background 

In this section we present the terminology and notation relevant to chemical 
reaction networks and the main results from the literature upon which we 
will be building. 

2.1 Chemical Reaction Networks 

We will consider the chemical species or reactants of a network to be given 
by the set S = {Xi, X2, . . . , Xn}. The combined elements on the left-hand 
and right-hand side of a reaction are given by linear combinations of these 
species. These combined terms are called complexes and will be denoted by 
the set C = {Ci, C2, . . . , Cm} where 



c^ = T. 



a-ijXj, i = 1, . . . ,m 



and the Oij are nonnegative integers called the stoichiometric coefficients. 
We define the reaction set to be 7^ = {(Q, Cj) \ Ci reacts to form Cj} where 
the property {Ci,Cj) € 7^ will more commonly be denoted Ci — )■ Cj. To 
each {Ci,Cj) £ TZ we will associate a positive rate constant kij > and to 
each {Ci,Cj) ^ IZ we will set kij = 0. The triplet M = {S,C,TZ) will be 
called the chemical reaction network. 

The above formulation naturally gives rise to a directed graph C{V, E) 
where the set of vertices is given by F = C, the set of directed edges is 
given hy E = TZ, and the rate constants kij corresponds to the weights of 
the edges from d to Cj . In the literature this has been termed the reaction 



graph of the network J20j. Since complexes may be involved in more than 
one reaction, as a product or a reactant, there is further graph theory we 
may consider. A linkage class is a maximally connected set of complexes, 
that is to say, two complexes are in the same linkage class if and only if there 
is a sequence of reactions in the reaction graph (of either direction) which 
connects them. A reaction network is called reversible if Ci — t- Cj for any 
Ci,Cj £ C implies Cj -^ Ci. A reaction network is called weakly reversible 
if Ci — )• Cj for any Ci,Cj £ C implies there is some sequence of complexes 
such that Ci = C^(i) -^ C^(2) -^ > ^i^ii-i) ^ C'mW = ^r 

A directed graph is called strongly connected if there exists a directed 
path from each vertex to every other vertex. A strongly connected component 
of a directed graph is any set of vertices for which paths exists from each 
vertex in the set to every other vertex in the set. For a weakly reversible 
network, the linkage classes clearly correspond to the strongly connected 
components of the reaction graph. 

Assuming mass-action kinetics, the dynamics of the specie concentra- 
tions over time is governed by the set of differential equations 

f = F-A.-*(x) (1) 

where x = [xi 0:2 • • • 2;„] is the vector of reactant concentrations. The 
stoichiometric matrix Y contains entries [Y]ij = Uji and the Kirchhoff or 
kinetics matrix A^ is given by 

[Akh = \ -^^1.'^^^^^' ^[ '=^. (2) 

I kji if ? / J. 

Finally, the mass-action vector ^(x) is given by 

n 

^j{^) = llxP\ j = l,...,m. (3) 

4 = 1 

2.2 Sparse and Dense Realizations 

Under the assumption of mass-action kinetics, it is possible for two reaction 
networks J\f^^' and A/"'^^ to give rise to the same set of governing differential 
equations. In other words, it is possible that 



yd) . 41) . .i/(i)(x) = y(2) . ^(2) . ^(2)(^) ^ f(^)^ V X G M^ 



>0 



where Y^^', A^ , and ^(*)(x), i = 1, 2, are the stoichiometric matrices, kinet- 
ics matrices, and mass-action vectors defined for M^^' and M^'^' respectively. 



The networks M^^' and M^'^' are called different realizations of the kinetics 
f (x) although it will sometimes be more convenient to consider J\f^^' as an 
alternative realization of M^"^' or vice- versa. 

In [30] the author presents an algorithm for producing sparse and dense 
alternative realizations of a given network M' , i.e. realizations with the 
fewest and greatest number of reactions capable of generating the same 
kinetics ([l]) as that given by M' . Key to the analysis is fixing the matrix 
Y to contain only the (source or product) complexes corresponding to the 
network M' . The problem of finding an alternative realization A^ of A/"' then 
becomes one of finding a kinetics matrix A^ such that 

y-Afc-M/(x) = y-^'fc-vi/(x). 

If we set M = Y ■ A'l^ and impose that A^ be a kinetics matrix, dynamical 
equivalence can be guaranteed by the conditions 



(DE) 



Y-Ak = M 



m 



Y,[Mij = o, j = i, 



i=l 



(4) 



[Ak]ij > 0, i,j = l,...,m, i^ j 
[Aklii < 0, i 



A sparse (respectively, dense) realization is given by a matrix A/^ satis- 
fying Q with the most (respectively, least) off-diagonal entries which are 
zeroes. A correspondence between the non-zero off-diagonal entries in A^ 
and a positive integer value can be made by considering the binary variables 
Sij G {0, 1} which will keep track of whether a reaction is 'on' or 'off', i.e. 
we have 

6ij = 1 ^ [Aklij > e, i, j = 1, . . . , m, i / j 

for some sufficient small < e ^ 1, where the symbol '-f-T-' denotes the logical 
relation 'if and only if. These proposition logic constraints for the structure 
of a network can then be formulated as the following linear mixed-integer 
constraints (see, for example, [25j ) : 

{0 < [Ak]ij - e6ij, i,j = l,...,m, i y^ j 
0< -[Ak]ij +Uij6ij, i,j = l,...,m, i^j (5) 

5ije {0,1}, i,j = l,...,m,i^ j, 

where Uij > for i,j = 1, . . . ,m,i ^ j. The number of reactions present 
in the network corresponding to A^ is then given by the sum of the 5jj's so 



that the problem of determining a sparse network corresponds to satisfying 
the objective function 

m 

(Sparse) < minimize > 5ij (6) 

over the constraint sets Q and ([5| . Finding a dense network corresponds to 
maximizing the same function, which can also be stated as a minimization 
problem as 



(Dense) < 



minimize 2, ~^ij- C^) 



2.3 Weakly Reversible Networks 

Weakly reversible networks are a particularly important class of reaction 
networks because strong properties are known about their dynamics. Under 
a supplemental condition, which is easily derived from the reaction graph 
alone, it is known that there is a unique positive equilibrium concentra- 
tion within each invariant space of the network and that that equilibrium 



concentration is at least locally asymptotically stable (15 18 201. 



In 33 the authors introduce an algorithm for determining dense weakly 
reversible realizations of a given kinetics. The algorithm is based on the fact 
that there are no cycles involving elements in different strongly connected 
components of a reaction network [61, and that for a fixed complex set the 
structure of the dense realization of a network is unique and contains the 
structures of all other possible realizations as sub-graphs [32]. Omitting 
technical details, the algorithm can be summarized as: 

1. Set the matrices Y and M and initialize /C = {}. 

2. Compute a dense realization A^ forcing the edges in /C to be excluded. 

3. Check whether Af. is weakly reversible (if so, end algorithm and return 
Ak). 

4. Find all edges in A^ which lead from one strongly connected compo- 
nent to another and add them to /C. 

5. Check whether these edges may be removed (if so, repeat steps (2)-(4); 
if not, end algorithm and return A/^ = 0). 



The algorithm has the drawbacks that it can only compute dense real- 
izations and not sparse ones, and that it requires potentially multiple MILP 
optimizations which are known to be NP-hard. In Section [3?T] we will present 
a method for determining both dense and sparse weakly reversible realiza- 
tions in a single MILP optimization step. 

3 Original Results 

In this section, we extend the results of [33] by showing how the requirement 
of weak reversibility can be formulated as a linear constraint. Consequently, 
the question of determining a sparse or dense weakly reversible realization 
can now be handled in a single MILP optimization step. We also extend 
this framework to cover the linearly conjugate networks introduced in [27| . 

3.1 Weak Reversibility as a Linear Constraint 

In this section we show that the requirement of weak reversibility can be 
formulated as a linear constraint. We require the following classical result 
about weakly reversible networks, which is modified from Theorem 3.1 of (16] 
and Proposition 4.1 of |14|: 



Theorem 3.1. Let A^ be a kinetics matrix and let Aj, i = 1, . . . ,i, denote 
the support of the i^^ linkage class. Then the reaction graph correspond- 
ing to Ak is weakly reversible if and only if there is a basis of ker{Ak), 
{b(i),...,bW}, such that, fori = !,...,£, 



An immediate consequence of Theorem 3.1 is that there is a vector b G 



M!!^Q n ker(Afc) if and only if the reaction graph corresponding to A^ is 
weakly reversible. In other words, we can guarantee weak reversibility by 
imposing the condition 

Akh = (8) 

for some b G M™q. This is a nonlinear constraint in the /cjj's and 6j's. In 
order to make it linear, we consider the matrix A/^ with entries 

[Ak]ij = [Ak]ij • bj. (9) 



It is clear from ([9j) that Aj^ encodes a kinetics matrix and that 1 £ M™ (the 
m-dimensional vector containing only ones) lies in ker(ylfc). Moreover, it is 



easy to see that A^ encodes a weakly reversible network if and only if A^ 
corresponds to a weakly reversible network. We can therefore check weak 
reversibility of the chemical reaction network corresponding to Ak with the 
linear conditions 



(WR') < 



k]ij=0, j = '^,---,m 



i=l 
m 

[Aklij > 0, i,j = l,...,m, iy^ j 



m 



(10) 



[Ak]ii<0, i = l,...,m. 

By solving for the diagonal elements of Ak, the set of constraints (10) can 
be simplified to 

(WR) { 



[Aklij > 0, ij = l,...,m, iy^j. 



,m 



(11) 



No condition comparable to y • ^fc = M exists for the matrix A^ so that 
we are left to optimization with respect to the internal entries of both Ak 
and Ak- Given appropriate choices of < e <C 1 and Uij > 0, i,j = 1, . . . ,m, 
i ^ j, we can impose 



(WR-S) 



0<[Ak 



eSi, 



l,...,m, ij^j 



^k\ij ^Oij , I, J — i,...,..D, " 7- J ('\'}\ 

< -[Ak]ij + Uij5ij, i,j = l,...,m, i^ j 



as well as (l5j) to ensure that both Ak and Ak contain zero and non-zero 
entries in the same places so that they correspond to reaction graphs with 
the same structure. 



3.2 Linear Conjugacy 



In [21] the authors extended the concept of dynamical equivalence to linear 
conjugacy. In their framework, two networks M and M' are said to be 
linearly conjugate if there is a linear mapping which takes the flow of one 
network to the other. The case of two networks being realizations of the 
same kinetics is encompassed as a special case of linear conjugacy taking 
the transformation to be the identity. (For a more complete introduction to 
the notion of dynamical equivalence and conjugacy see [23 or |34j.) 



Importantly, linearly conjugate networks share the same qualitative dy- 
namics (e.g. number and stability of equilibria, persistence/extinction of 
species, dimensions of invariant spaces, etc.). Similarly with different re- 
alizations of the same kinetics (II]), if a network with unknown kinetics is 
linearly conjugate to a network with known dynamics, then the qualitative 
properties of the second network are transferred to the first. 



The main result of 21 is the following. We have adopted the notation 
to match that contained in this paper. The notation is sufficiently distinct 
that we will prove the result independently here. 

Theorem 3.2. Consider two mass-action systems M = {S,C,TZ) and M' = 
{S,C',TZ') and let Y be the stoichiometric m,atrix corresponding to the com- 
plexes in either network. Consider a kinetics matrix A^ corresponding to N 
and suppose that there is a kinetics matrix A}, with the same structure as N' 
and a vector c S M"g such that 

Y-Ak = T-Y-At (13) 

where T =diag{c}. Then J\f is linearly conjugate to J\f' with kinetics matrix 

A'k = Ak-dtag{^{c)}. (14) 

Proof. Let $(xo,t) correspond to the flow of (IT]) associated to the reaction 
network A/". Consider the linear mapping h(x) = T~^ • x where T =diag{c}. 
Now define ^(yo, t) = T'^ ■ $(xo, t) so that $(xo, t) = T ■ $(yo, t). 
Since $(xo,t) is a solution of (jTj) we have 

l>'(yo,t) = T-i-cD'(xo,t) 

= T-i-y-^fc-^($(xo,t)) 

= T-i-r-y-A-^(r-l>(yo,t)) 

= y-A-diag{^(c)}-^(l>(yo,t)). 



It is clear that $(yo,t) is the flow of M corresponding to the reaction net' 



work A/"' with the kinetics matrix given by (14). We have that h($(xo, t)) = 
<&(h(xo),t) for all xq G M"q and t > where yo = h(xo) since yo = 
|)(yo, 0) = T^i • $(xo, 0) = T-i • XQ. It fohows that the networks M and J\f' 
are linearly conjugate and we are done. D 



While the results of 21 give conditions for two networks to be linearly 



conjugate, and therefore exhibit the same qualitative dynamics, no general 



methodology is provided for determining linearly conjugate networks when 
only a single network is provided. 

In cases where the dynamics of a network J\f is suspected to behave 
like a weakly reversible network, it is beneficial to extend the optimization 



algorithm outlined in Section 3.1 to linearly conjugate networks. This can 
be accomplished by replacing the set of constraints Ml) with 



(LC) 



Y ■Ab = T-^ 


■M 








)_,[Abh = 0, 


J -- 


= 1, 


. . . ,m 




[Ab]ij > 0, i 
[Ab]ii<0, i 
e<Cj< 1/e, 


J = 

= 1, 

j = 


1,- 
-- 1, 


m 
..,n 


/j 



(15) 



where M = Y ■ A/., T =diag{c}, and < e ^ 1, and replacing the set of 
constraints ([s]) by 



< [Ab]ij - eSij, i,j = l,...,m, i y^ j 
(LC-S) { < -[Ab]ij + Uijdij, i,j = l,...,m, i^ j 

6ijG{0,l}, i,j = l,...,m,i^j, 



(16) 



where Uij > for i,j = l,...,m,i^j. 

Ab has the same structure as the kinetics matrix A'f^ corresponding to the 
conjugate network, and this matrix has the same structure as the matrix A^ 
given by ^ (replacing A^ by A'^^). Consequently, the problem of determining 
a sparse or dense weakly reversible network which is linearly conjugate to 
a given kinetics can be given by optimizing either (rol) or ([7]), respectively, 
over the constraint sets (15), (16), (11), and (12). The kinetics matrix A'/^ 



for the linearly conjugate network is given by (14). 



4 Examples 

In this section we will consider two examples from the literature which 



demonstrate how the MILP optimization algorithm outlined in Section 3.2 



is capable of efficiently finding sparse and dense weakly reversible networks 
which are linearly conjugate to a given network M. We also consider one new 
example which illustrates how the algorithm is capable of finding networks 
with linearly conjugate dynamics for which no trivial linear conjugacy exists. 



10 



Example 1: Consider the chemical reaction network M given by 

Xi + 2X2 ^ Xi 

2X1 + X2 — > SX2 

Xi + 3X2 — > Xi + X2 — y 3X1 + X2 ■ 



M 



This network was first considered in 21 where the authors showed that 



it was hnearly conjugate to a specified weakly reversible network for all 

values of a > 0. It was further analysed with the value a = 1.5 in [33] where 

the authors found a dense weakly reversible realization through successive 

MILP optimizations. We will reproduce this result using our one-step MILP 

algorithm and also produce a sparse realization. 

We have 

"112 113 

2 13 3 11 



Y 



and 



M 






-2 


2 


-3 


2 


-2 



and set e = 1/a = 2/3 and Uij = 20, i,j = 1,...,7, i 7^ j. The MILP 
problem for a dense weakly reversible linearly conjugate network, possibly 
accounting for a non-trivial linear conjugacy mapping, is 





minimize 


i,3 = l 


over the constraint set 






Y-Ai, 


= T-^.M 




m 
1=1 


= 0, j = i,.. 


. ,m 


m 


m 




)l. [^^]^^- 


= > ; [AkU 


3 = 


«=i,«7^i 


1=1, i^j 





JlJ 



,m 



< [Ab\ij -e-6ij,i,j = l,...,7,iy^ j 
< -[Ablij + Uij ■dij,i,j = l,...,7,i^ j 
< [Aklij -e-6ij,i,j = l,...,7,iy^ j 
< -[Ak]ij + Uij ■6ij,i,j = l,...,7,i^ j 



11 



where T =diag{c}, and the decision variables 

[Ah]ij > 0, [Ak]ij > 0, for i,j = l,...,7,i^ j 
[Ab]ii<0, fori = l,...,7 
£ < Ci < 1/e, for i = 1, 2 

6ij G {0, 1} , for i, j = 1, . . . , 7, i / j. 



Solving for Ai, with GLPK and applying ( 14 ) gives the kinetics matrix 



A, 
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and values ci = 1,C2 = 1 (i.e. the linear transformation is the identity). 
The network structure is given graphically in Figure ffla). Although the 
rate constants differ due to differing bounds, this has the same network 



structure as the dense weakly reversible network obtained in 33 



A sparse weakly reversible network is generated by optimizing 



minimize 






over the same constraint set. Solving for A}, with GLPK with the bound 
e 



0.1 and applying (14) gives the kinetic matrix 



Ak 
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-100 








10 
































100 





-500 








50 














-10 



























and values ci = 10 and C2 = 5. This is therefore an example of a network 
with a non-trivial linear conjugacy and corresponds to the weakly reversible 
network given in Figure p[b). 



12 



11/3 



(a) X1+2X2 ^^ X1+X2 (b) X1+2X2 ^^ X1+X2 



2/3 



:2/3 



2/3 



t^ 



2/3IJ2 



^/3 

X^+3X2 < 2/3 2X^4X2 



'^ 



500 



10 



X,+3X2 ^^ 2X,+X2 



Figure 1: Dense (a) and sparse (b) weakly reversible networks which are 
linearly conjugate to J\f. 

Example 2: Consider the kinetics scheme 



Xl = X3 — X1X2 + X3X4 — 2x1X2X3 
X2 = X^ — X1X2 + 2x3X4 — 4x1X2X3 

X3 = —2x3 + X1X2 — X1X2X3 + 2x4 
X4 = X1X2 — X3X4 + 4x1X2X3 — 3x4 



first considered in [31]. Using the algorithm given in |17 and 31 
determine a kinetic realization involving the complexes 

Ci = 2X3, C2 = X3 + X4, C3 = Xi + 2X3, C4 = X2 + 2X3, 

C5 = X-i, Cq = Xi + Xs + X4, Cj = X2 + X3 + X4 

Cs = Xi + X2,Cg = Xi + 2X2 + X3, do = Xi, Cii = X2 

C\2 = Xi + X2 + X4, Cl3 = Xi + X2 + X3, Cl4 

Ci5 = X\ + 2X2, C16 = Xi + 2X2 + X3 + X4 
Ci7 = 3X4, C18 = X3 + 3X4, C\ 



(17) 



we can 



: 2X2 + X3 



^9 



S-3 
2X4. 



In 133 the authors use the algorithm given in Section 2.3 to determine 



a dense weakly reversible realization. The algorithm required three MILP 
optimizations, three searches for strongly connected components, and took 
80.5s to complete. Carrying out either the MILP optimization algorithm 



outlined in Section 3.2 for a dense or for a sparse weakly reversible network, 



with bounds e = 0.1 and Uij = 10, i,j = l,..., 19, we arrive at the solution 

kl8 = k29 = k82 = 0.1, kg2 = kg(n) = 0.001, A:(i7)i = 0.01, 
Ci = C2 = C3 = C4 = 0.1 



13 



and the rest of the entries zero (the transformation is a scahng of the iden- 
tity). This corresponds to the network given in Figure |4] which has the same 
network structure as the networks obtained in [31] and [33]. Our algorithm 
was able to obtain the answer in a single MILP optimization step and took 
less than a tenth of a second to compute. (The difference in rate constants 
occurs as a result of the scaling of concentration variables permitted by lin- 
ear conjugacy.) 



0.01 



2X. 



0.1 



> X1+X2 



,0.1 



0.001 



3X4 < Xi+2X2+X3^ X3+X4 

^ 0.001 ^ ^ "* ^Ta 



Figure 2: Weakly reversible realization of the kinetics (17). This realization 
is both dense and sparse. 



Example 3: Consider the kinetics scheme 



Xl 
X3 



2^13^2 



2Xx + 2:1X3 



2 2, 2 

X ^ 0Ct2 \ Xl Xo 



(18) 



™2 
Xl 



3X1X3. 



Using the algorithm given in llTj and 131], we can determine a kinetic real- 
ization involving the complexes 



Ci 



Xl + 2X2, C2 = 2X1 + 2X2, C3 = 2X1 + X2, 
2X1 ,€5 = Xi,Ce = 2X1 + X3 , C7 = Xl + 2X3 
2X1 + 2X3, Cg = Xl + X2 + 2X3, Cio = Xl + X3. 



With this fixed complex set, we can carry out the MILP optimization 
procedure outlined in Section [3. 2| to find sparse and dense weakly reversible 



networks which are linearly conjugate to a network with kinetics (18). We 

have 

"12 2 2 12 12 11 

2210000010 

0000012221 



Y 



14 



and 



M 



1 








-2 


10 





-1 








10 











1 


-3000 



With the bounds e = 1/20 and Uij = 20 for i,j = 1,...,10, i ^ j, 
the algorithm gives us the sparse network given in Figure Bfa) (conjugacy 
constants ci = 20, C2 = 2, and cs = 5) and the dense network given in 
Figure Qb) (conjugacy constants ci = 20/3, C2 = 20/33, and C3 = 5/3). 
It is interesting to note that the sparse and dense networks utihze different 
complexes and that the ratio of conjugacy constants differ between the sparse 
and dense networks. It is worth noting that the sparse realization is also 
deficiency zero so that the Deficiency Zero Theorem can be applied [15,18 
20 . Consequently, solutions of ( 18 ) satisfy all of the stringent dynamical 



restrictions typically reserved for complex balanced systems. 



(a) X1+2X2 

125\ 



25 



2X1+2X2 

^400 



X1+2X3 <=^ 2X1 

' -^ ^ 40 ' 



(b) X1+2X2 °-^ 2X1+2X2 




- 0.926s ^,, 

X,+2X^^=± 2X1 



2X1 +X2 



Figure 3: Weakly reversible networks which are linearly conjugate to a net- 



work with the kinetics (18). The network in (a) is sparse while the network 
in (b) is dense. The parameter values in (b) have been rounded to three 
significant figures. 



5 Conclusions 



In this paper we have proposed an algorithm for determining linearly con- 
jugate weakly reversible realizations of reaction networks. In contrast to 
the method presented in (33| , the present approach is based on the well- 
known fact that the kernel of the Kirchhoff matrix of weakly reversible 
networks always contains a strictly positive vector. The main advantages of 



our algorithm compared to 33 are the following. Firstly, linear conjugacy 



15 



theory [21] has been included into the optimization framework, and the pa- 
rameters of the corresponding hnear coordinates transformation belong to 
the set of unknowns. Secondly, our algorithm requires only one MILP step, 



and therefore it is numerically significantly more effective than 33 if the 
problem dimension is similar to what is shown in the examples. Thirdly, ad- 
ditional structural constraints such as density or sparsity of the solution net- 
work can be directly included into the optimization problem. The presented 
results clearly contribute to further widening the application possibilities of 
the known strong results in chemical reaction network theory. 

There are still several very important questions which remain to be an- 
swered: 

1. While the algorithm is effective and efficient for finding alternative net- 
works for a given kinetics, we are often interested in questions which 
can be answered for all kinetics satisfying certain initial structural 
properties (i.e. general rather than specied rate constants). The algo- 
rithm is currently unable to answer such questions. 

2. Many dynamical properties are known for systems satisfying network 
structure properties not included in weak reversibility theory (8l[9l[T2] . 
Extending the optimization framework to include these results would 
greatly expand the scope of chemical reaction networks with known 
dynamics. 
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